import geopandas as gpd
import pandas as pd  # 新增：解决pd.isna()未定义问题
import trimesh
import numpy as np
import os
from shapely.ops import unary_union
from shapely.validation import make_valid, explain_validity
from shapely import is_valid
from shapely.errors import TopologicalError, GEOSException

def qifengshan_farmland_to_3d(geojson_path, output_path, scale=1000000, export_format='glb', extrusion_height=1.0):
    """
    将七峰山林场耕地GeoJSON（CRS84/EPSG:4326）转换为Unity可用的3D模型
    核心命名规则：小地块名称 = LIN_CHANGM（林场名） + LIN_BAN2（林班） + XIAO_BAN2（小班）
    例：LIN_CHANGM=七峰山林场、LIN_BAN2=1078、XIAO_BAN2=00121 → 命名为「七峰山林场107800121」

    参数说明：
    - geojson_path: 输入七峰山林场耕地GeoJSON文件路径
    - output_path: 输出3D模型路径（默认GLB格式，Unity兼容性最佳）
    - scale: 坐标缩放系数（经纬度放大倍数，默认1e6，避免Unity中坐标过小）
    - export_format: 导出格式（支持glb/gltf/obj/stl/ply，默认glb）
    - extrusion_height: 3D模型挤压高度（默认1.0，根据Unity场景尺度调整）
    """
    # ---------------------- 1. 基础输入校验 ----------------------
    # 检查输入文件存在性
    if not os.path.exists(geojson_path):
        print(f"【错误】输入GeoJSON文件不存在：{geojson_path}")
        return

    # 创建输出目录（不存在则新建）
    output_dir = os.path.dirname(output_path)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
        print(f"【提示】已创建输出目录：{output_dir}")

    # 验证导出格式合法性
    supported_formats = ['glb', 'gltf', 'obj', 'stl', 'ply']
    if export_format not in supported_formats:
        print(f"【错误】不支持的导出格式：{export_format}，自动切换为默认格式glb")
        export_format = 'glb'

    # 修正输出路径扩展名（确保格式匹配）
    if not output_path.lower().endswith(f'.{export_format}'):
        output_path = f"{output_path}.{export_format}"
        print(f"【提示】调整输出路径为：{output_path}")

    # 碰撞检测库可用性检查（python-fcl，用于精确防重叠）
    try:
        import fcl
        fcl_available = True
        print("【提示】检测到python-fcl库，将使用精确碰撞检测（避免网格重叠）")
    except ImportError:
        fcl_available = False
        print("【提示】未检测到python-fcl库，使用轻量化边界框碰撞检测")

    # ---------------------- 2. 读取并解析七峰山林场耕地GeoJSON ----------------------
    try:
        # 读取CRS84（EPSG:4326）坐标系统的GeoJSON数据（保留所有字段）
        gdf = gpd.read_file(geojson_path, crs="EPSG:4326")
        print(f"\n【成功】读取七峰山林场耕地GeoJSON：共包含 {len(gdf)} 个地块，坐标系统→CRS84（EPSG:4326）")
        print(f"【字段信息】数据包含 {len(gdf.columns)} 个属性字段，将完整保留为网格元数据")
    except Exception as e:
        print(f"【错误】读取GeoJSON失败：{str(e)}（请确认文件为标准FeatureCollection格式）")
        return

    # 检查核心命名字段（必须包含：LIN_CHANGM、LIN_BAN2、XIAO_BAN2；唯一标识：MZGUID）
    core_naming_fields = ["LIN_CHANGM", "LIN_BAN2", "XIAO_BAN2", "MZGUID"]
    missing_fields = [field for field in core_naming_fields if field not in gdf.columns]
    if missing_fields:
        print(f"【错误】GeoJSON缺少核心命名/标识字段：{missing_fields}，无法按规则生成小地块名称")
        return

    # 处理字段空值（避免命名异常，空值替换为"未知"）
    for field in core_naming_fields:
        gdf[field] = gdf[field].fillna("未知")
        # 确保数值型字段转为字符串（如LIN_BAN2=1078→"1078"，避免拼接时丢失格式）
        if gdf[field].dtype in [np.int64, np.float64]:
            gdf[field] = gdf[field].astype(str).str.strip()
        else:
            gdf[field] = gdf[field].astype(str).str.strip()

    # ---------------------- 3. 修复几何图形拓扑错误 ----------------------
    print("\n【处理】开始修复几何图形拓扑错误（自相交、空洞、重复顶点等）...")
    valid_geoms = []
    for idx, geom in enumerate(gdf.geometry):
        # 跳过空几何或无效几何
        if geom is None or geom.is_empty:
            valid_geoms.append(None)
            continue

        try:
            # 检查几何有效性，无效则修复
            if not is_valid(geom):
                # 打印前10个无效几何的原因（避免日志冗余，仅显示关键信息）
                if idx < 10:
                    err_info = explain_validity(geom)[:50] + "..." if len(explain_validity(geom))>50 else explain_validity(geom)
                    mzguid = gdf.iloc[idx]['MZGUID'][:8] + "..."  # 隐藏部分GUID，保护数据隐私
                    print(f"【修复】地块{idx}（MZGUID:{mzguid}）拓扑错误：{err_info}")
                # 修复无效几何（支持MultiPolygon/Polygon，优先保留面积最大子区）
                fixed_geom = make_valid(geom)
                if fixed_geom.geom_type == "MultiPolygon":
                    fixed_geom = max(fixed_geom.geoms, key=lambda p: p.area if p.exterior else 0)
                valid_geoms.append(fixed_geom)
            else:
                valid_geoms.append(geom)
        except Exception as e:
            mzguid = gdf.iloc[idx]['MZGUID'][:8] + "..."
            print(f"【错误】地块{idx}（MZGUID:{mzguid}）几何修复失败：{str(e)}")
            valid_geoms.append(None)

    # 更新有效几何并过滤空值
    gdf['geometry'] = valid_geoms
    gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty]
    if len(gdf) == 0:
        print("【错误】所有地块几何修复后均无效，无法生成3D模型")
        return
    print(f"【成功】几何修复完成：保留 {len(gdf)} 个有效耕地地块")

    # 检查原始数据拓扑重叠（提示数据质量，避免Unity中显示重叠）
    try:
        union_geom = unary_union(gdf.geometry)
        if union_geom.geom_type in ["MultiPolygon", "GeometryCollection"]:
            print(f"【提示】检测到 {len(union_geom.geoms)} 个独立区域，原始数据可能存在重叠或不连通")
    except (TopologicalError, GEOSException) as e:
        print(f"【警告】拓扑重叠检查失败：{str(e)}，跳过该检查")

    # ---------------------- 4. 坐标预处理（经纬度→平面偏移） ----------------------
    # 收集所有有效顶点（仅保留X/Y坐标，移除可能的Z值，确保Unity中平面定位准确）
    all_coords = []
    valid_geometry_indices = []  # 记录有效地块在gdf中的索引
    for idx, row in gdf.iterrows():
        geom = row['geometry']
        # 统一处理MultiPolygon和Polygon类型（适配数据几何格式）
        if geom.geom_type == "MultiPolygon":
            polygons = list(geom.geoms)
        elif geom.geom_type == "Polygon":
            polygons = [geom]
        else:
            mzguid = row['MZGUID'][:8] + "..."
            print(f"【警告】地块（MZGUID:{mzguid}）几何类型为{geom.geom_type}，仅支持Polygon/MultiPolygon，已跳过")
            continue

        # 提取多边形顶点（闭合多边形需至少4个点：起点=终点，中间至少2个点，避免生成无效网格）
        poly_coords = []
        for poly in polygons:
            if poly.exterior is None:
                continue
            coords = list(poly.exterior.coords)
            coords_2d = [coord[:2] for coord in coords]  # 强制转为2D坐标（移除Z值）
            if len(coords_2d) >= 4:
                poly_coords.extend(coords_2d)

        # 记录有效坐标和地块索引
        if poly_coords:
            all_coords.extend(poly_coords)
            valid_geometry_indices.append(idx)
        else:
            mzguid = row['MZGUID'][:8] + "..."
            print(f"【警告】地块（MZGUID:{mzguid}）有效顶点不足4个，已跳过")

    # 检查是否收集到有效坐标（无有效坐标则终止）
    if not all_coords:
        print("【错误】未收集到任何有效顶点数据，无法生成3D模型")
        return

    # 计算坐标偏移量（以最小经纬度为原点，避免Unity中坐标值过大导致精度问题）
    all_coords_np = np.array(all_coords, dtype=np.float64)
    min_lon, min_lat = all_coords_np.min(axis=0)  # CRS84：X=经度，Y=纬度
    offset = np.array([min_lon, min_lat], dtype=np.float64)
    print(f"\n【坐标配置】偏移原点（经度：{min_lon:.8f}，纬度：{min_lat:.8f}），缩放系数：{scale}")

    # ---------------------- 5. 构建3D网格（核心：专属命名规则） ----------------------
    scene = trimesh.Scene()  # Trimesh场景（整合所有3D网格，便于统一导出）
    valid_meshes = []        # 存储所有有效3D网格
    existing_meshes = []     # 碰撞检测缓存（避免网格重叠，确保Unity显示正常）

    # 配置三角化引擎（提升多边形转网格的效率和精度，优先使用第三方引擎）
    try:
        import mapbox_earcut
        os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "earcut"
        print("【提示】使用高效三角化引擎：mapbox_earcut（处理复杂多边形更快）")
    except ImportError:
        try:
            import triangle
            os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "triangle"
            print("【提示】使用三角化引擎：triangle")
        except ImportError:
            print("【提示】未找到第三方三角化引擎，使用Trimesh默认引擎（适合简单多边形）")

    # 核心参数配置（可根据数据复杂度和Unity性能需求调整）
    SIMPLIFY_TOLERANCE = scale * 0.00001  # 多边形简化阈值（减少顶点数，降低Unity渲染压力）
    MIN_VERTICES = 3                      # 简化后最小顶点数（低于则保留原始数据，避免网格失效）
    MIN_DISTANCE_THRESHOLD = 1e-5         # 碰撞检测最小距离（微调重叠网格的位置）

    # 遍历有效地块生成3D网格（核心：按LIN_CHANGM+LIN_BAN2+XIAO_BAN2命名）
    for idx in valid_geometry_indices:
        row = gdf.iloc[idx]
        # 提取核心命名字段（生成小地块名称）
        lin_changm = row["LIN_CHANGM"].strip()  # 林场名（如：七峰山林场）
        lin_ban2 = row["LIN_BAN2"].strip()      # 林班（如：1078）
        xiao_ban2 = row["XIAO_BAN2"].strip()    # 小班（如：00121）
        mzguid = row["MZGUID"].strip()          # 唯一标识（用于重名兜底）

        # 核心：小地块命名规则 → LIN_CHANGM + LIN_BAN2 + XIAO_BAN2（如：七峰山林场107800121）
        base_mesh_name = f"{lin_changm}{lin_ban2}{xiao_ban2}"
        # 重名兜底：若不同地块命名重复，补充MZGUID后8位（避免Unity中网格名称冲突）
        if len([m for m in valid_meshes if m.name.startswith(base_mesh_name)]) > 0:
            base_mesh_name = f"{base_mesh_name}_{mzguid[-8:]}"

        # 重新获取当前地块的多边形（处理MultiPolygon包含多个子区的情况）
        geom = row["geometry"]
        polygons = list(geom.geoms) if geom.geom_type == "MultiPolygon" else [geom]

        # 处理每个子多边形（MultiPolygon的每个子区单独生成网格，避免合并失真）
        for sub_idx, poly in enumerate(polygons):
            # 子多边形命名（多子区时添加序号，如：七峰山林场107800121_SUB0）
            mesh_name = f"{base_mesh_name}_SUB{sub_idx}" if len(polygons) > 1 else base_mesh_name

            # 检查子多边形外边界有效性（无外边界则跳过）
            if poly.exterior is None:
                print(f"【警告】{mesh_name} 无外边界，已跳过")
                continue
            coords = list(poly.exterior.coords)
            coords_2d = [coord[:2] for coord in coords]  # 强制2D坐标（移除Z值）

            # 步骤1：坐标转换（经纬度→平面坐标，应用偏移和缩放，适配Unity尺度）
            planar_coords = (np.array(coords_2d, dtype=np.float64) - offset) * scale

            # 步骤2：简化复杂多边形（顶点数>500时执行，平衡精度和性能）
            path = trimesh.load_path(planar_coords)
            original_vertex_count = len(planar_coords)
            if original_vertex_count > 500:
                simplified_path = path.simplify(tolerance=SIMPLIFY_TOLERANCE)
                simplified_count = len(simplified_path.vertices)
                if simplified_count >= MIN_VERTICES:
                    path = simplified_path
                    print(f"【简化】{mesh_name}：顶点数 {original_vertex_count} → {simplified_count}（提升Unity渲染效率）")
                else:
                    print(f"【提示】{mesh_name} 简化后顶点不足{MIN_VERTICES}个，保留原始顶点")

            # 检查简化后顶点数（避免生成无效网格）
            if len(path.vertices) < MIN_VERTICES:
                print(f"【警告】{mesh_name} 有效顶点不足{MIN_VERTICES}个，已跳过")
                continue

            # 步骤3：生成3D网格（挤压2D多边形为立体模型，Z轴方向为高度）
            try:
                # 挤压生成3D网格（height=挤压高度，即模型厚度）
                mesh = path.extrude(height=extrusion_height)
                # 处理Trimesh部分版本挤压返回列表的特殊情况
                if isinstance(mesh, list):
                    mesh = mesh[0] if mesh else None
            except Exception as e:
                print(f"【警告】{mesh_name} 挤压生成网格失败，尝试手动构建：{str(e)}")
                # 手动构建平面网格（备选方案，避免完全失败）
                vertices = np.column_stack([planar_coords, np.zeros(len(planar_coords))])  # Z轴设为0（基础平面）
                faces = [[0, i, i+1] for i in range(1, len(planar_coords)-1)]  # 扇形三角化（确保网格闭合）
                if not faces:
                    print(f"【警告】{mesh_name} 无法手动构建网格，已跳过")
                    continue
                mesh = trimesh.Trimesh(vertices=vertices, faces=np.array(faces))

            # 检查网格有效性（顶点数和面数均需大于0，否则跳过）
            if mesh is None or len(mesh.vertices) == 0 or len(mesh.faces) == 0:
                print(f"【警告】{mesh_name} 未生成有效网格，已跳过")
                continue

            # 步骤4：碰撞检测（避免网格重叠，确保Unity中每个地块独立显示）
            collision_found = False
            if existing_meshes:
                if fcl_available:
                    # 精确碰撞检测（基于python-fcl，检测网格实体是否重叠）
                    if check_collision_fcl(mesh, existing_meshes):
                        # 沿X/Y轴微调位置（微小距离，不影响整体布局）
                        mesh.apply_translation([MIN_DISTANCE_THRESHOLD, MIN_DISTANCE_THRESHOLD, 0])
                        print(f"【碰撞】{mesh_name} 与已有网格重叠，已微调位置（偏移量：{MIN_DISTANCE_THRESHOLD}）")
                        collision_found = True
                else:
                    # 轻量化边界框碰撞检测（速度快，适合大规模数据）
                    current_bbox = mesh.bounds  # 当前网格边界框（[min坐标, max坐标]）
                    for existing_mesh in existing_meshes:
                        existing_bbox = existing_mesh.bounds
                        # 轴对齐边界框碰撞检测（判断X/Y轴方向是否重叠）
                        if (current_bbox[0][0] < existing_bbox[1][0] and
                                current_bbox[1][0] > existing_bbox[0][0] and
                                current_bbox[0][1] < existing_bbox[1][1] and
                                current_bbox[1][1] > existing_bbox[0][1]):
                            mesh.apply_translation([MIN_DISTANCE_THRESHOLD, MIN_DISTANCE_THRESHOLD, 0])
                            print(f"【碰撞】{mesh_name} 与已有网格边界框重叠，已微调位置")
                            collision_found = True
                            break

            # 步骤5：设置网格元数据（完整保留GeoJSON所有字段，便于Unity中查询耕地属性）
            mesh_metadata = {"mesh_name": mesh_name, "CRS": "EPSG:4326",
                             "offset_longitude": round(min_lon, 8), "offset_latitude": round(min_lat, 8),
                             "scale_factor": scale, "extrusion_height": extrusion_height}
            # 遍历GeoJSON所有字段，添加到元数据（确保属性不丢失）
            for col in gdf.columns:
                if col != "geometry":  # 跳过几何字段（已用于生成网格）
                    val = row[col]
                    # 处理特殊值（空值→"未知"，数值→保留2位小数，避免科学计数法）
                    if pd.isna(val):
                        mesh_metadata[col] = "未知"
                    elif isinstance(val, (int, float)):
                        mesh_metadata[col] = round(val, 2) if val is not None else "未知"
                    else:
                        mesh_metadata[col] = str(val).strip() if val is not None else "未知"
            mesh.metadata = mesh_metadata
            mesh.name = mesh_name  # Unity编辑器中显示的网格名称（便于识别）

            # 步骤6：添加到场景和缓存（用于后续碰撞检测和统一导出）
            valid_meshes.append(mesh)
            existing_meshes.append(mesh)
            print(f"【生成】{mesh_name}：顶点数{len(mesh.vertices)}，面数{len(mesh.faces)}（元数据包含{len(mesh_metadata)}个字段）")

    # ---------------------- 6. 导出3D模型（Unity兼容格式） ----------------------
    if not valid_meshes:
        print("\n【错误】未生成任何有效3D网格，模型导出失败")
        return

    try:
        # 导出3D模型（包含所有网格和完整元数据）
        scene.add_geometry(valid_meshes)
        scene.export(file_type=export_format, file_obj=output_path)
        print(f"\n【🎉 导出成功】七峰山林场耕地3D模型生成完成！")
        print(f"文件路径：{output_path}")
        print(f"核心统计：{len(valid_meshes)}个3D网格，对应{len(valid_geometry_indices)}个有效耕地地块")
        print(f"命名规则：小地块名称 = LIN_CHANGM + LIN_BAN2 + XIAO_BAN2（如：七峰山林场107800121）")
        print(f"格式说明：{export_format}（Unity导入推荐GLB格式，单个文件包含所有网格+元数据，无需额外依赖）")
    except Exception as e:
        print(f"\n【错误】{export_format}格式导出失败：{str(e)}")
        # 备选格式兜底（优先尝试OBJ格式，兼容性最强，适合所有3D软件）
        fallback_format = "obj" if export_format != "obj" else "glb"
        fallback_path = output_path.rsplit(".", 1)[0] + f".{fallback_format}"
        try:
            scene.export(file_type=fallback_format, file_obj=fallback_path)
            print(f"【备选导出成功】路径：{fallback_path}（格式：{fallback_format}，可正常导入Unity）")
        except Exception as fe:
            print(f"【错误】备选格式{fallback_format}导出也失败：{str(fe)}")


def check_collision_fcl(mesh, existing_meshes):
    """
    基于python-fcl的精确碰撞检测
    判断当前网格是否与已有网格集合中的任意网格重叠（避免Unity中显示重叠）
    返回：True（重叠）/False（不重叠）
    """
    try:
        import fcl
        # 构建当前网格的BVH（层次包围盒）碰撞模型（高效检测）
        current_bvh = fcl.BVHModel()
        current_bvh.beginModel(len(mesh.vertices), len(mesh.faces))
        current_bvh.addSubModel(mesh.vertices, mesh.faces)
        current_bvh.endModel()
        current_obj = fcl.CollisionObject(current_bvh)

        # 遍历已有网格，逐一检测碰撞
        for existing_mesh in existing_meshes:
            existing_bvh = fcl.BVHModel()
            existing_bvh.beginModel(len(existing_mesh.vertices), len(existing_mesh.faces))
            existing_bvh.addSubModel(existing_mesh.vertices, existing_mesh.faces)
            existing_bvh.endModel()
            existing_obj = fcl.CollisionObject(existing_bvh)

            # 执行碰撞检测（获取碰撞结果）
            collision_req = fcl.CollisionRequest()
            collision_res = fcl.CollisionResult()
            fcl.collide(current_obj, existing_obj, collision_req, collision_res)

            # 若检测到碰撞，立即返回True
            if collision_res.is_collision:
                return True
        # 所有已有网格均无碰撞，返回False
        return False
    except Exception as e:
        print(f"【警告】FCL碰撞检测出错：{str(e)}，默认判定为无碰撞（可忽略，不影响模型生成）")
        return False


# ---------------------- 执行入口（修改路径后直接运行） ----------------------
if __name__ == "__main__":
    # 1. 配置输入输出路径（请根据您的文件实际位置修改！）
    INPUT_GEOJSON = r"C:\Users\ma\Desktop\山河屯数据\耕地区域\七峰山林场耕地.geojson"  # 输入GeoJSON路径
    OUTPUT_MODEL = r"C:\Users\ma\Desktop\山河屯数据\耕地区域\七峰山林场耕地_3D模型.glb"  # 输出3D模型路径

    # 2. 调用转换函数（参数可根据Unity场景需求调整）
    qifengshan_farmland_to_3d(
        geojson_path=INPUT_GEOJSON,
        output_path=OUTPUT_MODEL,
        scale=1000000,          # 坐标缩放系数（经纬度×1e6，适合Unity显示，无需修改）
        export_format='glb',    # 导出格式（Unity推荐glb，如需OBJ可改为'obj'）
        extrusion_height=1.2    # 3D模型高度（如1.2单位，可根据场景尺度调整，如改为2.0）
    )
